Attention select de mass (dépendance de mixOmics) prend le pas sur dplyr
tCPMG_vector <- unlist(as.vector(t(CPMG)))
CPMG_vector <- unlist(as.vector(CPMG))
tNOESY_vector <- unlist(as.vector(t(NOESY)))
NOESY_vector <- unlist(as.vector(NOESY))
par(mfrow = c(2, 2))
hist(CPMG_vector, breaks = 100, xlab = 'CPMG', main = "Histogram of CPMG" )
hist(log2(CPMG_vector), breaks = 100, xlab = 'log2 CPMG', main = "Histogram of log2 CPMG")
hist(NOESY_vector, breaks = 100, xlab = 'Noesy', main = "Histogram of NOESY")
hist(log2(NOESY_vector), breaks = 100, xlab = 'log2 Noesy', main = "Histogram of log2 NOESY")df_row_sum <- apply(NOESY,MARGIN = 1, FUN= mean)
df_col_sum <- apply(NOESY,MARGIN = 2, FUN= mean)
meanr <- mean(df_row_sum)
meanc <- mean(df_col_sum)
minr <- min(df_row_sum)
maxr <- max(df_row_sum)
minc <- min(df_col_sum)
maxc <- max(df_col_sum)
medr <- median(df_row_sum)
medc <- median(df_col_sum)
df <- data.frame(
nb_ech = c(846,215),
moy = c(meanr,meanc),
median = c(medr,medc),
min = c(minr,minc),
max = c(maxr, maxc)
)
rownames(df) <- c("Echantillon","Bucket")
kable(df, align = "c", caption = "Tableau de moyenne des échantillons et buckets pour NOESY") %>%
kable_styling(
bootstrap_options = "striped",
full_width = F,
position = "center")| nb_ech | moy | median | min | max | |
|---|---|---|---|---|---|
| Echantillon | 846 | 0.0066202 | 0.0065891 | 0.0059730 | 0.0090406 |
| Bucket | 215 | 0.0066202 | 0.0026412 | 0.0000208 | 0.2255073 |
standardisation par la médiane et mise à l’échelle par la plage inter-quartile
tNOESY <- t(NOESY)
log2_tNOESY <- log2(tNOESY)
log2_NOESY <- log2(NOESY)
###Standardisation par la médiane
stdr_medians <- apply(log2_tNOESY, 2, median)
series_median <- median(stdr_medians)
log2_NOESY_centered <- data.frame(matrix(
nrow = nrow(log2_tNOESY),
ncol = ncol(log2_tNOESY)))
colnames(log2_NOESY_centered) <- colnames(log2_tNOESY)
rownames(log2_NOESY_centered) <- rownames(log2_tNOESY)
for (j in 1:ncol(log2_tNOESY)) {
log2_NOESY_centered[, j] <- log2_tNOESY[, j] - stdr_medians[j] + series_median
}
### Mise à l'échelle
sampleIQR <- apply(log2_NOESY_centered, 2, IQR)
seriesIQR <- median(sampleIQR)
scalingFactors <- sampleIQR / seriesIQR
NOESY_standard <- data.frame(matrix(
nrow = nrow(log2_tNOESY),
ncol = ncol(log2_tNOESY)))
colnames(NOESY_standard) <- colnames(log2_NOESY_centered)
rownames(NOESY_standard) <- rownames(log2_NOESY_centered)
for (j in 1:ncol(log2_NOESY_centered)) {
NOESY_standard[, j] <-
(log2_tNOESY[, j] - stdr_medians[j] ) / scalingFactors[j] + series_median
}
par.ori <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
par(mar = c(5, 5, 2, 1))
#par(mar = c(5, 6, 4, 1)) ou 4 6 5 1
sample_size <- 30
## select sample indices
selected_samples <- sort(sample(
x = 1:ncol(tNOESY),
size = sample_size,
replace = FALSE))
boxplot(tNOESY[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "original values",
xlab = "value",
ylab = "sample")
boxplot(log2_tNOESY[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "log2-transformed",
xlab = "log2(value)",
ylab = "sample")
boxplot(log2_NOESY_centered[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "Median-based centered",
xlab = "log2(value)",
ylab = "sample")
boxplot(NOESY_standard[,selected_samples],
horizontal = TRUE,
las = 1,
yaxt="n",
main = "Standardized\n(median centering, IQR scaling)",
xlab = "log2(value)",
ylab = "sample")Les 2 premiers axes expliquent plus de 46% de la variance exprimée.
NOESY_stdr <- t(NOESY_standard)
res.pca <- PCA(NOESY_stdr, graph = FALSE)
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 40))Quelque soit le groupe étudié (sexe, âge, tabac), sur les axes 1 et 2, les différents points se superposent, il n’est donc pas possible de distinguer un groupe d’un autre.
L’individu L36R338 semble excentré sur l’axe3 (axe1-20) par rapport aux autres points et dans une moindre mesure V351R356 qui lui ne ressort pas sur les 2 premiers axes.
Contribution des 20 buckets qui contribuent le plus à la variance expliquée.
Cette contribution est >1
fviz_pca_biplot(res.pca, repel = TRUE, axes = c(1,2),
select.var = list(contrib = 20),
geom.ind = "point",
fill.ind = Sample_ID$SEX,
pointshape = 21, pointsize = 2,
palette = "jco",
addaEllipses = TRUE,
col.var = "contrib",
gradient.cols = "BuGn",
legend.title = list(fill = "Sexe", color = "Contrib"))Juste pour montrer le code mais les résultats sont identiques. Remarque : Chaque axe est symétrique donc qu’un individu soit à-20ici et +20 avec factominer revient au même
trans.pca <- pca(NOESY_stdr, ncomp=5, center = TRUE, scale = TRUE)
plotIndiv(trans.pca, comp = c(1,2),ind.names = FALSE,
group = as.numeric(as.factor(Sample_ID$SEX)),
legend.title = "sexe")Regarde si on peut observer naturellement une clusterisation en fonction d’une catégorie
Quelque soit la catégorie sexe, age tabac, nous n’observons aucune clusterisation.
res.spca.NOESY <- spca(NOESY_stdr, scale=TRUE,
ncomp=3, keepX=c(10,10,10))
plotIndiv(res.spca.NOESY, ind.names = FALSE,
group=Sample_ID$SEX,
pch = as.numeric(Sample_ID$SEX),
pch.levels =Sample_ID$SEX,
legend = TRUE)Clusturisation par PLS-DA structure de correlation entre CPMG and NOESY variables séléctionné sur la composante 1
Cluster_CPMG_NOESY
trans.plsda <- plsda(NOESY_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotIndiv(trans.plsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur NOESY")Ne sont représentés uniquement les buckets qui contribuent à au moins 50%
plotVar(trans.plsda, var.names = FALSE, cutoff = 0.5,
title = "Correlation circle plot avec cutoff 0.5")Possibilité de faire de la prédiction dans le cas ou les ind seraient en sous groupes, ce qui n’est pas notre cas.
ceci est juste un test mais pourrait être interressant afin de prédire une variable à partir de 100aines que nous avons
trans.testplsda <- plsda(NOESY_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
background <- background.predict(trans.testplsda, comp.predicted=2,
dist = "max.dist")
plotIndiv(trans.testplsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur NOESY", background = background,
legend.title="Sexe")Les variables sélectionnées pour la composante 1 sont exprimées de manière positives pour le sexe1 et négative pour le sexe2.
Pour l’âge : la composante est positive pour les tranches d’âge de 50 à 69et négative pour les 20-29 ans.
Pour le statut tabagique la composante est exclusivement négatives pour le statut non et ex fumeur et positif pour les fumeurs.
MyResult.plsda2 <- plsda(NOESY_stdr, Sample_ID$SEX, ncomp=5)
plotLoadings(MyResult.plsda2, contrib = 'max', method = 'mean')MyResult.plsda2 <- plsda(NOESY_stdr, Sample_ID$age, ncomp=5)
plotLoadings(MyResult.plsda2, contrib = 'max', method = 'mean')Chez les fumeurs la contribution des variables a un effet positif alors qu’elle est négative dans les 2 autres catégories.
res.plsda <- plsda(NOESY_stdr, Sample_ID$tabac, ncomp=5)
plotLoadings(res.plsda, contrib = 'max', method = 'mean', comp=1)Résultats simlaires à la PLSDA
trans.splsda <- splsda(NOESY_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE, legend=TRUE, legend.title="sexe",
ellipse = TRUE, title="Sparse PLS-DA sur NOESY")trans.splsda <- splsda(NOESY_stdr, Sample_ID$AGE ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE, legend=TRUE, legend.title="age",
ellipse = TRUE, title="Sparse PLS-DA sur NOESY")trans.splsda <- splsda(NOESY_stdr, Sample_ID$tabac ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE,legend=TRUE, legend.title="Tabac",
ellipse = TRUE, title="Sparse PLS-DA sur NOESY")Représentation par sexe uniquement
Ne sont représentés uniquement les buckets qui contribuent à au moins 50%
trans.splsda <- splsda(NOESY_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotVar(trans.splsda, var.names = FALSE, cutoff = 0.5,
title = "Correlation circle plot avec cutoff 0.5")Représentation des 10 premières contributions
MyResult.splsda2 <- splsda(NOESY_stdr, Sample_ID$SEX, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Sexe")MyResult.splsda2 <- splsda(NOESY_stdr, Sample_ID$tabac, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Tabac")MyResult.splsda2 <- splsda(NOESY_stdr, Sample_ID$AGE, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Age")